#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern double cvt_density();
extern double smooth_front();
extern void handle_triple();

extern double laplace_above,laplace_below,laplace_middle,
              density_above,density_below,density_middle;

/* compute the viscous forces
   input ix,iy
   output cfl_laplace - time restriction due to explicit central diff. approx 
     of the laplacian.
   ut,vt - values are UPDATED (already contain values from convective forces)
     central differencing is used.
*/ 
void do_laplacian(cfl_laplace)
double *cfl_laplace;
{
double ux[2][2],uy[2][2],vx[2][2],vy[2][2],ur[2][2],
       stress1[2][2],stress2[2][2],stress3[2][2],
       laplacectr[2][2],d1ctr[2][2],d2ctr[2][2],d3ctr[2][2];
double density,rplus,rminus,r_val;
double laplace_center,viscous_u,viscous_v;
int m,n;
int ix,iy;

  *cfl_laplace = 0.0;
  if (laplace_coeff != 0.0)
   {
    for_all_points
     if (! point_ok)
      {
       u[ix][iy]=gu(ix,iy);
       v[ix][iy]=gv(ix,iy);
       d1[ix][iy]=gd1(ix,iy);
       d2[ix][iy]=gd2(ix,iy);
       d3[ix][iy]=gd3(ix,iy);
      }

    for_domain
    {
    
    r_val=x[ix];
    rplus=r_val+h/2.0;
    rminus=r_val-h/2.0;

#define dcenter(gdtype,ix1,iy1) ((gdtype[ix-1+ix1][iy-1+iy1]+ \
 gdtype[ix+ix1][iy-1+iy1]+gdtype[ix-1+ix1][iy+iy1]+gdtype[ix+ix1][iy+iy1])/4.0)

    for (m=0;m<2;m++)
     for (n=0;n<2;n++)
      {
       ux[m][n]=(u[ix+m][iy+n]+u[ix+m][iy+n-1]-u[ix+m-1][iy+n]-
                 u[ix+m-1][iy+n-1])/(2.0*h);
       vx[m][n]=(v[ix+m][iy+n]+v[ix+m][iy+n-1]-v[ix+m-1][iy+n]-
                 v[ix+m-1][iy+n-1])/(2.0*h);
       uy[m][n]=(u[ix+m][iy+n]-u[ix+m][iy+n-1]+u[ix+m-1][iy+n]-
                 u[ix+m-1][iy+n-1])/(2.0*h);
       vy[m][n]=(v[ix+m][iy+n]-v[ix+m][iy+n-1]+v[ix+m-1][iy+n]-
                 v[ix+m-1][iy+n-1])/(2.0*h);
       

       if (use_axisymmetric)
       {
       if ((m==0)&&(fabs(rminus)<almost_same))
        ur[m][n]=ux[m][n];
       else
        {
         ur[m][n]=(u[ix+m][iy+n]+u[ix+m][iy+n-1]+u[ix+m-1][iy+n]+
                   u[ix+m-1][iy+n-1])/4.0;
         if (m==0) ur[m][n] /= rminus;
         else ur[m][n] /= rplus;
        }
       }  /* only compute for axisymmetric r-z flow */


       d1ctr[m][n]=dcenter(d1,m,n);
       d2ctr[m][n]=dcenter(d2,m,n);
       d3ctr[m][n]=dcenter(d3,m,n);

       if (use_bmo)
        handle_triple(&(d1ctr[m][n]),&(d2ctr[m][n]),&(d3ctr[m][n]),
           viscosity_spread*h);

       laplacectr[m][n]=smooth_front(d1ctr[m][n],d2ctr[m][n],d3ctr[m][n],
         laplace_above,laplace_below,laplace_middle,viscosity_spread); 

       stress1[m][n]=2.0*laplacectr[m][n]*ux[m][n];
       stress2[m][n]=laplacectr[m][n]*(vx[m][n]+uy[m][n]);
       stress3[m][n]=2.0*laplacectr[m][n]*vy[m][n];
      }

    density=cvt_density(d1[ix][iy],d2[ix][iy],d3[ix][iy]);
    laplace_center=smooth_front(d1[ix][iy],d2[ix][iy],d3[ix][iy],
       laplace_above,laplace_below,laplace_middle,viscosity_spread);

    if (use_axisymmetric)
     {
      /*
      viscous_u=(laplace_coeff/density)*(
       (rplus*(stress1[1][1]+stress1[1][0])-
        rminus*(stress1[0][1]+stress1[0][0]))/(2*r_val*h)-
       2.0*laplace_center*u[ix][iy]/sqr(r_val)+
       (stress2[1][1]+stress2[0][1]-stress2[1][0]-stress2[0][0])/(2.0*h)
      );
      */
      viscous_u=(laplace_coeff/density)*(
       (stress1[1][1]+stress1[1][0]-stress1[0][1]-stress1[0][0])/(2*h)+
       2.0*laplace_center*(ur[1][1]+ur[1][0]-ur[0][1]-ur[0][0])/(2*h)+
       (stress2[1][1]+stress2[0][1]-stress2[1][0]-stress2[0][0])/(2*h)
      );
       

      viscous_v=(laplace_coeff/density)*(
       (rplus*(stress2[1][1]+stress2[1][0])-
        rminus*(stress2[0][1]+stress2[0][0]))/(2*r_val*h)+
       (stress3[1][1]+stress3[0][1]-stress3[1][0]-stress3[0][0])/(2.0*h)
      );
     }
    else
     {
      viscous_u=(laplace_coeff/density)*(
       (stress1[1][1]+stress1[1][0]-stress1[0][1]-stress1[0][0])/(2.0*h)+
       (stress2[1][1]+stress2[0][1]-stress2[1][0]-stress2[0][0])/(2.0*h)
      );

      viscous_v=(laplace_coeff/density)*(
       (stress2[1][1]+stress2[1][0]-stress2[0][1]-stress2[0][0])/(2.0*h)+
       (stress3[1][1]+stress3[0][1]-stress3[1][0]-stress3[0][0])/(2.0*h)
      );
     }

    ut[ix][iy]+=viscous_u;
    vt[ix][iy]+=viscous_v;

    if (fluid_select==222)
     printf("ix,iy,vu,vv%5d%5d%15.10f%15.10f\n",ix,iy,viscous_u,viscous_v);
    }   /* for_domain */
    
    *cfl_laplace=fabs( (14.0/3.0)*(laplace_coeff/density)*laplace_center/h);
   }   /* laplace_coeff != 0.0 */

}   /* do_laplacian */
